# 11/9 file based on 221014; robustness check of implications for not constraining race options to ancestry

#Set Directory
setwd('/Users/drschaef/Google Drive/Projects/Adriana/Raw data/Longitudinal Data/')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Raw data/Longitudinal Data/')

library(network)
library(sna)
library(statnet)
library(dplyr)
library(RSiena)

load("~/Google Drive/Projects/Adriana/Papers/Homophily/allDatPreSIENA221014.RData")

# function to create heritage homophily matrices
makeHeritageSim <- function(ATT) {
	rDums <- as.matrix(ATT[,which(names(ATT)=='A'): which(names(ATT)=='N')], ncol=5)
	nInCommon <- rDums %*% t(rDums)
	N <- dim(ATT)[1]
	nRace <- rowSums(rDums)
	maxInCommon <- pmax(matrix(rep(nRace,N),ncol=N), t(matrix(rep(nRace,N),ncol=N)) )
	maxInCommon[maxInCommon==0] <- .000001
	homophFull <- nInCommon / maxInCommon
	homophPart <- homophFull
	homophPart[homophPart==1] <- 0
	homophPart[homophPart>0] <- 1
	homophFull[homophFull<1] <- 0
	diag(homophFull) <- 0
	diag(homophPart) <- 0
	# separate into multi-multi dyads and multi-mono dyads
	momoFull <- homophFull * outer(ATT$mono==1, ATT$mono==1, '*')
	mumuFull <- homophFull * outer(ATT$multi==1, ATT$multi==1, '*')
	mumuPart <- homophPart * outer(ATT$multi==1, ATT$multi==1, '*')
	mumoPart <- homophPart * (outer(ATT$mono==1, ATT$multi==1, '*') + 
							 outer(ATT$multi==1, ATT$mono==1, '*'))
	mumu <- (homophFull + homophPart) * outer(ATT$multi==1, ATT$multi==1, '*')
	mumo <- homophPart * outer(ATT$multi==1, ATT$mono==1, '*')
	momu <- homophPart * outer(ATT$mono==1, ATT$multi==1, '*')
	return(array(c(homophFull, homophPart, momoFull, mumuFull, mumuPart, mumoPart, mumu, mumo, momu), dim=c(N,N,9)))
	}

# a function to create the SIENA object
# as part of the function, create heritage homophily dyadic covariate (scaled using .6)
#  - 	HOM <- HOM1[,,1] + .6*HOM1[,,2] # uses hSim layers 1 and 2

makeSienaData <- function (NET1, NET2, ATT, ECA, RACE1, RACE2, POPSIZE, STARTWAVE, DVS=NULL) {
	# calculate homophily on heritage
	hSim <- makeHeritageSim(ATT)
	erSimScaled <- hSim[,,1] + .6*hSim[,,2]
	erSimNomomo <- hSim[,,4] + hSim[,,5] + hSim[,,6]

	# counts
	nActors <- dim(NET1)[1]
	nRace <- dim(RACE1)[2]
	
	# define different node sets:
	students <- sienaNodeSet(nActors, nodeSetName="students")
	erGroups <- sienaNodeSet(nRace, nodeSetName="erGroups")

	# dependent networks
	friendship <- sienaNet(array( c(NET1, NET2), dim = c(nActors, nActors, 2)), nodeSet="students")
	raceNet <- sienaNet(array( c(RACE1, RACE2), dim = c(nActors, nRace, 2)), "bipartite", nodeSet=c("students","erGroups"))

	# endogenous version of ECA (rigged to be exogenous)
	ecaEndog <- sienaNet(array( c(ECA[[1]], ECA[[2]]), dim = c(nActors, nActors, 2)), nodeSet="students")

	# student attributes
	female <- coCovar(ATT$Female, nodeSet="students")
	gradecohort <- coCovar(ATT$Gradecohort, nodeSet="students")
	pared <- coCovar(ATT$ParEd, nodeSet="students")
	pared8 <- coCovar(ATT$pared8, nodeSet="students")
	immgen <- coCovar(ATT$ImmigGen, nodeSet="students")
	multi <- coCovar(ATT$multi, center=F, nodeSet="students")
	mono <- coCovar(ATT$mono, center=F, nodeSet="students")
	nis <- coCovar(ATT$NIS, nodeSet="students")
	perla <- coCovar(ATT$PERLA, nodeSet="students")
	racedif <- coCovar(ATT$RACEDIF, nodeSet="students")
	egoA <- coCovar(ATT$A, nodeSet="students")
	egoB <- coCovar(ATT$B, nodeSet="students")
	egoL <- coCovar(ATT$L, nodeSet="students")
	egoW <- coCovar(ATT$W, nodeSet="students")
	egoN <- coCovar(ATT$N, nodeSet="students")
	egoO <- coCovar(ATT$O, nodeSet="students")
	egoM <- coCovar(ATT$MUL, nodeSet="students")
	
	# exogenously changing covariates (treated as constant within multigroup setup)
	# - academic performance
	tempAcad <- ATT[,grepl("academic", names(ATT), fixed=T)]
	tempAcad[tempAcad==-9] <- NA
	tempAcad[is.na(tempAcad[,1]),1] <- tempAcad[is.na(tempAcad[,1]),2]  # if t1 missing, impute with t2
	acad <- coCovar(tempAcad[,STARTWAVE], nodeSet="students")
	tempAcad9 <- ATT[,grepl("acad9", names(ATT), fixed=T)]
	tempAcad9[is.na(tempAcad9[,1]),1] <- tempAcad9[is.na(tempAcad9[,1]),2]  # if t1 missing, impute with t2
	acad9 <- coCovar(tempAcad9[,STARTWAVE], nodeSet="students")
	# - exogeneous version of current race
	tempRace <- data.frame(ATT$chooseOne1, ATT$chooseOne2, ATT$chooseOne3)
	raceExog <- coCovar(tempRace[,STARTWAVE], nodeSet="students")
		
	# - discrimination, overall
	temp5 <- data.frame(ATT$disc1, ATT$disc2, ATT$disc3)
	disc <- coCovar(temp5[,STARTWAVE], nodeSet="students")
	
	# race attributes: other = 6 (reference category)
	asian <- coCovar(c(1,0,0,0,0,0,0), center=F, nodeSet='erGroups')
	black <- coCovar(c(0,1,0,0,0,0,0), center=F, nodeSet='erGroups')
	latin <- coCovar(c(0,0,1,0,0,0,0), center=F, nodeSet='erGroups')
	white <- coCovar(c(0,0,0,1,0,0,0), center=F, nodeSet='erGroups')
	natam <- coCovar(c(0,0,0,0,1,0,0), center=F, nodeSet='erGroups')
	orace <- coCovar(c(0,0,0,0,0,1,0), center=F, nodeSet='erGroups')
	mrace <- coCovar(c(0,0,0,0,0,0,1), center=F, nodeSet='erGroups')
	pPop <- coCovar(POPSIZE, nodeSet='erGroups')

	# dyadic covariates (student X student)
	ecaDum <- coDyadCovar(ECA[[1]], center=F, nodeSets=c("students", "students"))
	herRaceSim <- coDyadCovar(erSimScaled, center=F, nodeSets=c("students", "students"))
	herRaceSim2 <- coDyadCovar(erSimNomomo, center=F, nodeSets=c("students", "students"))
	mumu <- coDyadCovar(hSim[,,7], center=F, nodeSets=c("students", "students"))
	mumo <- coDyadCovar(hSim[,,8], center=F, nodeSets=c("students", "students"))
	momu <- coDyadCovar(hSim[,,9], center=F, nodeSets=c("students", "students"))

	# dyadic covariates (student X race)
	immgenA <- coDyadCovar(outer(immgen, asian, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenB <- coDyadCovar(outer(immgen, black, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenL <- coDyadCovar(outer(immgen, latin, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenW <- coDyadCovar(outer(immgen, white, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenN <- coDyadCovar(outer(immgen, natam, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenO <- coDyadCovar(outer(immgen, orace, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenM <- coDyadCovar(outer(immgen, mrace, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaA <- coDyadCovar(outer(perla, asian, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaB <- coDyadCovar(outer(perla, black, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaL <- coDyadCovar(outer(perla, latin, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaW <- coDyadCovar(outer(perla, white, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaN <- coDyadCovar(outer(perla, natam, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaO <- coDyadCovar(outer(perla, orace, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaM <- coDyadCovar(outer(perla, mrace, "*"), center=T, nodeSets=c("students","erGroups"))
	discA <- coDyadCovar(outer(temp5[,STARTWAVE], asian, "*"), center=T, nodeSets=c("students","erGroups"))
	discB <- coDyadCovar(outer(temp5[,STARTWAVE], black, "*"), center=T, nodeSets=c("students","erGroups"))
	discL <- coDyadCovar(outer(temp5[,STARTWAVE], latin, "*"), center=T, nodeSets=c("students","erGroups"))
	discW <- coDyadCovar(outer(temp5[,STARTWAVE], white, "*"), center=T, nodeSets=c("students","erGroups"))
	discN <- coDyadCovar(outer(temp5[,STARTWAVE], natam, "*"), center=T, nodeSets=c("students","erGroups"))
	discO <- coDyadCovar(outer(temp5[,STARTWAVE], orace, "*"), center=T, nodeSets=c("students","erGroups"))
	discM <- coDyadCovar(outer(temp5[,STARTWAVE], mrace, "*"), center=T, nodeSets=c("students","erGroups"))

	myRaceFriendsData <- sienaDataCreate( friendship, raceNet,
	       female, gradecohort, pared, pared8, immgen, acad, acad9, multi, mono, nis, perla, racedif,
	       asian, black, latin, white, natam, orace, mrace, pPop,
	       disc, egoA, egoB, egoL, egoW, egoN, egoO, egoM,
	       ecaDum, herRaceSim, herRaceSim2, mumu, mumo, momu,  
	       nodeSets=list(students, erGroups))
	       
	return(myRaceFriendsData)
	
	}


# remove structural zeros from race matrices
raceM1_[raceM1_==10] <- 0
raceM2_[raceM2_==10] <- 0
raceM_2[raceM_2==10] <- 0
raceM_3[raceM_3==10] <- 0
raceS1_[raceS1_==10] <- 0
raceS2_[raceS2_==10] <- 0
raceS_2[raceS_2==10] <- 0
raceS_3[raceS_3==10] <- 0

rfnetM12 <- makeSienaData(matM1_, matM_2, attM12, ecaM12list, raceM1_, raceM_2, pPopM, 1, 'race+friendship')
rfnetM23 <- makeSienaData(matM2_, matM_3, attM23, ecaM23list, raceM2_, raceM_3, pPopM, 2, 'race+friendship')
rfnetS12 <- makeSienaData(matS1_, matS_2, attS12, ecaS12list, raceS1_, raceS_2, pPopS, 1, 'race+friendship')
rfnetS23 <- makeSienaData(matS2_, matS_3, attS23, ecaS23list, raceS2_, raceS_3, pPopS, 2, 'race+friendship')
rfnetM123_remove10s <- sienaGroupCreate(list(rfnetM12, rfnetM23))
rfnetS123_remove10s <- sienaGroupCreate(list(rfnetS12, rfnetS23))

# remove ALL structural values (zeros and ones)
raceM1_[raceM1_==11] <- 1
raceM2_[raceM2_==11] <- 1
raceM_2[raceM_2==11] <- 1
raceM_3[raceM_3==11] <- 1
raceS1_[raceS1_==11] <- 1
raceS2_[raceS2_==11] <- 1
raceS_2[raceS_2==11] <- 1
raceS_3[raceS_3==11] <- 1

rfnetM12 <- makeSienaData(matM1_, matM_2, attM12, ecaM12list, raceM1_, raceM_2, pPopM, 1, 'race+friendship')
rfnetM23 <- makeSienaData(matM2_, matM_3, attM23, ecaM23list, raceM2_, raceM_3, pPopM, 2, 'race+friendship')
rfnetS12 <- makeSienaData(matS1_, matS_2, attS12, ecaS12list, raceS1_, raceS_2, pPopS, 1, 'race+friendship')
rfnetS23 <- makeSienaData(matS2_, matS_3, attS23, ecaS23list, raceS2_, raceS_3, pPopS, 2, 'race+friendship')
rfnetM123_remove10s11s <- sienaGroupCreate(list(rfnetM12, rfnetM23))
rfnetS123_remove10s11s <- sienaGroupCreate(list(rfnetS12, rfnetS23))

setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')

saveRDS(rfnetM123_remove10s, file = 'rfnetM123_remove10s.rds')
saveRDS(rfnetS123_remove10s, file = 'rfnetS123_remove10s.rds')
saveRDS(rfnetM123_remove10s11s, file = 'rfnetM123_remove10s11s.rds')
saveRDS(rfnetS123_remove10s11s, file = 'rfnetS123_remove10s11s.rds')



